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1.0  INTRODUCTION 

Op  deal  transmission  through  ciouds  is  of  importance  in  optical  communications 
and  remote  sensing  applications.  In  both  dowr'mk  and  uplink  communications 
systems  the  transmission  and  reflection  of  signal  and  background  light  from  clouds 
is  of  importance  in  determining  the  strength  of  the  perceived  signal-to-noise  ratio. 
Of  equal  importance  is  the  spatial,  angular,  and  temporal  spreading  of  the  signal 
pulse  in  the  cloud.  That  is,  due  to  random  scattering  events  in  the  cloud  a 
collimated  signal  will  be  physically  wider,  decollimated,  _  ;d  arrive  at  the  receiver 
over  time  periods  exceeding  the  original  pulse.  All  of  these  can  impact  the  signal- 
to-noise  ratio.  Likewise,  in  remote  sensing,  a  thorough  understanding  of  the 
reflection  from  cloud  surface  is  required  in  order  to  ascertain  the  cloud  optical 
properties. 

This  report  is  part  of  an  ongoing  effort  at  NOSC  to  characterize  clouds.  In  the  past, 
most  of  the  effort  has  been  placed  on  the  characterization  of  stratus  clouds.  Ln  this 
connection,  analytic  models  have  been  developed  for  cloud  transmission  and 
simulations  have  been  carried  out  for  spatial  and  temporal  spreading.  Now 
attention  is  being  turned  to  finite  clouds.  The  transmission  and  spreading  (spatial, 
angular,  and  temporal)  are  all  of  interest  for  finite  clouds. 

A  simulation  was  developed  for  radiative  transfer  in  finite  clouds  in  order  to  gain 
insight  into  the  cloud  characteristics.  This  report  describes  the  finite  cloud  model 
and  its  verification  for  a  stratus  cloud.  The  results  are  compared  with  the  current 
deterministic  link  models. 

Section  2  contains  a  technical  description  of  the  model.  Section  3  contains  the 
results  of  cloud  simulations  and  comparison  with  other  models.  Conclusions  and 
recommendations  are  proffered  in  Section  4.  References  can  be  found  in  Section  5 

2.0  TECHNICAL  DISCUSSION 


2.1  Background 

The  simulation  of  radiative  transfer  in  clouds  is  carried  out  bv  the  Monte  Carlo 
method.  In  a  Monte  Carlo  computation  one  photon  at  a  time  is  followed  in  its 
three-dimensional  path  through  a  scattering  medium.  Its  fate  is  determined  bv 
suitable  probability  distributions  for  mean  free  path,  absorption,  scattering  angle, 
reflection  angle  from  a  surface,  absorption  at  a  surface  and  so  on.  The  photon  is 
followed  until  it  is  absorbed,  detected,  or  leaves  the  field  of  interest.  A  sufficient 
number  ot  such  photons  are  followed  until  a  picture  of  the  system  emerges.  A 
detailed  accounting  of  the  Monte  Carlo  method  is  given  in  Appendix  A. 
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Monte  Carlo  simulations  have  been  used  in  the  past  at  NOSC  to  develop  models  tor 
spatial  and  temporal  spreading  in  stratus  clouds.  These  models  consist  or  curve  tits 
to  simulation  data  and  are  necessarily  limited  in  applicability  to  the  ranges  ot  optical 
and  physical  properties  used  in  the  simulations.  The  curve  tits  currently  in  use  are 
not  built  upon  analytic  models  and  essentially  have  no  physical  basis.  Thev  have 
not  been  adequately  documented  and,  as  will  be  seen  later,  may  not  be  applicable  to 
communications  problems. 

For  the  treatment  of  finite  clouds,  a  “cubic”  (actually  parallelpiped)  cloud  model  has 
been  developed.  The  model  is  supplemented  by  the  inter^’ty  reference  method 
which  is  used  for  the  prediction  of  photon  detection  by  surface  receivers.  In  the 
following  sections  the  cubic  cloud  model  and  intensity  reference  method  are 
described.  The  application  of  these  models  to  the  prediction  of  signal,  background, 
cloud  transmission  and  optical  thickness,  and  Signal-to-noise  ratio  are  described. 

2.2  Cubic  cloud  model 

The  cubic  cloud  model  consists  of  a  parallelpiped  cloud  situated  above  the  ground. 
The  cloud  geometry  and  coordinate  system  are  shown  in  Figure  1.  The  three  sides 
of  the  parallelpiped  as  well  as  the  zenith  and  azimuth  of  the  incoming  photons  may¬ 
be  specified  independently.  The  model  is  then  completely  specified  by  the  optica! 
thickness  of  the  cloud  and  the  asymmetry  parameter.  The  model  uses  the  Henyev- 
Greenstein  phase  function  for  the  prediction  of  scattering  angle  and  Poisson 
statistics  for  the  mean  free  path  (see  Appendix  A  for  a  more  complete  description.) 

Photons  may  enter  the  cloud  in  from  one  to  three  faces,  depending  upon  the  source 
zenith  and  azimuth  For  each  face  that  receives  photons  a  complete  simulation  is 
performed.  Actually,  these  are  identical  calculations  with  the  appropriate 
coordinate  transformation  and  source  angle.  Moreover,  the  number  of  photons 
entering  each  face  is  adjusted  so  that  each  face  has  the  same  energy  flux. 

Thus,  the  model  appropriately  accounts  for  the  photons  entering  different  taces  and 
computes  the  fraction  of  photons  exiting  from  each  of  the  six  faces.  In  addition,  the 
mean  and  standard  deviation  of  the  exit  angle  cosine,  average  position  and  distance 
traveled  are  computed.  Table  1  shows  a  typical  set  of  input  and  output. 

In  order  to  isolate  the  cloud  behavior,  the  present  cubic  cloud  mode1  does  not 
consider  the  surrounding  atmosphere. 
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2.3  Intensity  Reference  Method 

The  cubic  cloud  model  can  be  applied  to  communications  problems  by  allowing  the 
collection  of  photons  at  a  receiver.  Typically,  the  receiver  is  sufficiently  far  from  the 
cloud  and  sufficiently  small  that  literally  millions  of  photons  must  be  generated  for 
each  one  reaching  the  receiver.  Clearly,  this  is  undesirable  because  of  the  large 
computation  time  required.  Instead,  recourse  is  made  to  an  adaptation  of  the 
intensity  reference  method  due  to  Meier,  Lee,  and  Anderson  (1978).  In  this 
technique,  the  probability  of  reaching  the  receiver  is  computed  and  stored  as  each 
photon  exits  the  cloud.  The  accumulated  probabilities  are  proportional  to  the 
intensity.  A  more  detailed  description  of  the  intensity  reference  method  is  given  in 
Appendix  B. 

The  intensity  reference  method  considers  the  probability  of  the  photon  reaching  a 
detector  as  a  consequence  of  a  scattering  event.  Figure  2  shows  a  schematic  of  the 

method.  Starting  at  the  scattering  event  the  photon  has  a  probability  of  P  ( Gs)dco  of 
being  scattered  in  the  direction  of  the  receiver  (if  it  is  in  the  receiver  field  of  view). 
Along  that  path  there  may  also  be  a  probability  of  the  photon  being  absorbed  or 
scattered  out  of  the  path.  In  addition,  if  the  photon  is  being  multiply-scattered  then 
a  weighting  function  would  also  be  applied  if  absorption  is  present. 

The  advantage  of  this  method  is  that,  in  contrast  to  the  low  likelihood  of  a  photon 
actually  hitting  receiver,  there  will  always  be  some  probability  of  it  actually 
occurring.  Thus,  summing  probabilities  greatly  improves  the  statistics,  since  each 
scattering  event  contributes  to  the  intensity.  In  addition,  the  probability  is  computed 
exactly  from  the  scattering  point  to  the  receiver;  no  other  approximations  are 
introduced.  Moreover,  the  technique  lends  itself  to  vignetted  receivers,  i.e., 
receivers  with  a  limited  field  of  view  (which  lowers  the  likelihood  of  a  photon 
hitting  the  receiver  even  further). 

In  the  present  adaptation  of  the  intensity  reference  method  the  probability  of 
reaching  the  receiver  is  computed  as  the  photon  exits  the  cloud.  Thus,  the  photon  is 
followed  through  the  cloud  until  it  exits,  at  which  point  the  probability  of  being 
scattered  toward  the  receiver  is  calculated  from  the  last  collision.  More  generally, 
the  probabilities  would  be  computed  as  the  photon  is  scattered  in  the  atmosphere. 

A  variation  of  the  intensity  reference  method  was  tried  assuming  a  Lambertian 
radiance  profile  at  the  cloud  base.  In  this  case,  the  probability  is  computed  from  the 
radiance  profile  rather  than  the  scattering.  The  results  so  obtained  were 
unsatisfactory  because  of  inconsistencies  with  the  globa:  .  loud  transmission.  This 
led  to  an  investigation  of  the  radiance  profile  beneath  a  cloud.  These  results  showed 
that  while  the  radiance  profile  beneath  is  a  thick  cloud  is  relatively  independent  of 
cloud  optical  thickness  and  source  zenith  it  is  not  Lambertian.  Rather,  it  is  skewed 
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to  favor  forward  scattered  photons.  These  results  are  discussed  in  Section  3.  At  any 
rate,  this  method  was  not  compatible  with  the  intensity  reference  method. 

The  intensity  reference  method  is  well  suited  to  this  problem  and  would  work 
equally  well  if  there  was  an  atmosphere.  The  success  of  the  method  is  due  largely  to 
the  far-field  approximation  being  satisfied.  Specifically,  the  receiver  is  small 
compared  with  the  range  and  the  solid  angle  may  be  considered  as  a  differential. 
When  this  assumption  breaks  down  the  phase  function  must  be  integrated  over  the 
solid  angle.  This  is  discussed  further  in  Appendix  B. 

2.4  Receivers 

Two  types  of  receivers  were  modeled  to  work  in  conjunction  with  the  cubic  cloud 
model.  The  first  is  a  flat  plate  receiver  with  a  cosine  response  function  over  a  90- 
degree  half-angle  field  of  view.  The  second  is  a  vignetted  receiver  with  a  cosine 
response  over  a  specified  field  of  view  (less  than  90  degrees).  These  are  referred  to  as 
"ideal"  and  "actual"  receivers,  respectively.  They  can  also  be  thought  of  as 
representing  radiometer  and  communications  receivers,  respectively.  The  results 
for  the  ideal  receiver  can  be  used  to  determine  the  cloud  optical  thickness  from  a 
simulation  just  as  measured  irradiance  is  used  from  a  radiometer.  The  results  for 
the  actual  receiver  can  be  used  to  determine  signal,  background,  and  signal-to-noise 
ratio.  These  applications  are  discussed  in  the  following  sections. 

2.5  Cloud  Transmission  Model 

The  method  whereby  cloud  optical  thickness  can  be  inferred  from  a  measured 
surface  irradiance  is  well  known  (see,  for  example,  Waldman,  1986).  The  optical 
thickness  is  expressed  implicitly  in  the  physical  relationship  between  the  perceived 
irradiance,  cloud  transmission,  and  source  irradiance.  In  the  simulation  this 
translates  to 


Ar 


(1) 


where  xb  is  the  number  of  background  photons  reaching  the  receiver  of  area  Ar 
from  a  source  of  N  photons  at  a  direction  cosine  ^  on  a  cloud  of  area  Ac .  The 

optical  thickness  x  can  be  obtained  from  Eq.  (1)  by  iteration.  Lc  is  the  cloud 
transmission  given  by  King  and  Harshvardhan  (1986)  and  Waldman  (1987).  In 
developing  Eq.  (1)  an  ideal  (flat  plate)  receiver  has  been  assumed.  The  method 
applies  equally  well  to  a  vignetted  receiver  if  Eq.  (1)  is  modified  as  follows 
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*>  * 
Ar(l-cos20)  "  At 

where  6  is  the  receiver  field  of  view  half-angle. 

2.6  Signal-to-Noise  Model.' 

The  signal- to-noise  ratio  is  defined  here  as 


S /N  =  --T*— 

«,  +  nb  (3) 

where  n,  and  nb  are  the  number  of  signal  of  and  background  photoelectrons, 
respectively. 

For  a  downlink  these  are  defined  as  follows 

"■  =  W-fa’  (4) 

Note  that  E,Arfsg  in  Eq.  (4)  is  simply  the  fraction  of  signal  energy  reaching  the 

receiver  and  PbAr  fn  in  Eq.  (5)  is  the  fraction  of  background  power  reaching  the 
receiver.  Equations  (4)  and  (5)  take  the  following  form  for  a  simulation 


where  x,  and  xb  are  the  signal  and  background  photons  collected  in  the 
simulation,  respectively.  Notice  that,  in  general,  these  come  from  different 
simulations  because  the  signal  and  background  may  have  different  zenith  angles 
and  spot  sizes  on  top  of  the  cloud  (the  background  always  covers  the  entire  cloud, 
the  signal  spot  may  be  any  size). 
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Sigrtal-to-noise  ratio  predictions  from  the  simulations  can  be  compared  with  those 
from  a  deterministic  model.  Such  a  model  was  developed  for  this  purpose.  It 
contains  both  signal  and  background  calculations  for  a  receiver  below  a  stratus 
cloud.  Spot  shape  is  calculated  so  the  results  are  sensitive  to  receiver  position. 

There  is  no  atmosphere  in  the  model. 

2.7  Pulsewidth  Model 

The  receiver  integration  time,  T  is  modeled  as  it  is  in  the  deterministic  link  modei. 

For  a  matched  filter  it  is  usually  taken  as 

T  =  5tm  -  2,tjrf8  ^ 

Eq.  (8)  assumes  a  modified  gamma  function  pulse;  this  may  noc  generally  be  true. 
Appendix  C  contains  a  more  detailed  descirption  of  the  modified  gamma  function. 

The  receiver  integration  time  in  this  analysis  was  measured  from  the  simulation 
pulse  histogram. 

3.  RESULTS  AND  DISCUSSION 

Several  kinds  of  results  were  obtained  in  the  course  of  testing  and  verifying  the 
cubic  cloud  model.  These  are  described  below  along  with  some  typical  results  and 
their  interpretation. 

i 

3.1  Simple  cloud  model 

The  cubic  cloud  model  was  tested  on  its  own  prior  to  applying  it  to  the  signal-to- 
noise  calculation.  The  model  should  demonstrate  symmetry  and  global  energy 
conservation,  and  should  give  the  same  transmission  as  the  deterministic  model 
when  the  cloud  length  and  width  are  much  greater  than  height.  The  simulation 
was  also  tested  to  ensure  that  the  results  for  a  diffuse  source  agree  with  theoretical 
predictions;  to  wit,  that  a  diffuse  source  produces  the  same  result  as  an  equivalent 
direction  cosine  of  fj.  =  2/3 .  The  model  does  indeed  satisfy  all  these  criteria.  Tables 
1-3  shows  the  results  for  series  of  calculations  in  which  the  cloud  is  progressively 
increased  in  size  while  its  thickness  is  held  constant.  These  tables  show  that  as  the 
cloud  area  increases  the  photons  exiting  at  the  base  approach  the  value  of  cloud 
transmission  given  by  the  deterministic  model. 

John  Yen  of  NOSC  has  extended  this  cloud  model  to  study  the  irradiance  patterns 
on  the  surface  beneath  a  finite  cloud.  His  results  will  be  reported  separately. 
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3.2  Radiance  Profiles 

A  simplified  version  of  the  cloud  model  was  developed  to  studv  the  radiance 
profiles  upon  reflection  and  transmission  from  a  stratus  cloud.  The  simplifications 
consist  of  removing  the  cloud  edges  (i.e.,  infinite  cloud)  and  using  a  point  source 
.rather  than  one  distributed  over  the  cloud  area). 

The  program  produces  histograms  of  the  transmission  and  reflection  radiance 
profiles  expressed  as  probability  density'  functions  versus  angle.  Figures  3-8  show 
comparisons  of  the  reflection  results  with  those  for  a  Lambertian  profile  for  two 

optical  thicknesses  (r  =  20  and  100)  and  three  source  zeniths  (£-0°,  60°  and  diffuse). 

The  results  tor  the  source  at  zenith  (£= 0 °),  i.e.,  Figs.  3  and  6  show  non-Lambertian 
radiance  which  favors  smaller  angles,  this  effect  being  more  pronounced  at  higher 
optical  thickness.  The  results  for  an  oblique  source  (±=  60°),  i.e..  Figs.  4  and  7,  also 
exhibit  a  non-Lambertian  radiance  but  here  the  larger  angles  are  favored.  Bear  in 
mind,  however,  that  these  curves  represent  azimuthal  averages  and  are  not  really 
representative  of  any  particular  directional  reflectance  function.  The  results  for  a 
diffuse  source  (assumed  to  be  Lambertian)  are  shown  in  Figs.  5  and  S.  For  the 

smaller  optical  thickness  (r  =  20)  the  result  is  similar  to  that  for  the  oblique  source. 

For  the  larger  optical  thickness  (v  ~  100)  the  result  is  Lambertian. 

Figures  9-14  show  comparisons  of  the  transmission  results  with  those  for  a 
Lambertian  profile  for  the  same  conditions  as  above.  These  results  are  fairlv 
uniform  for  all  conditions.  Specifically,  all  these  cases  demonstrate  that  the  radiance 
profile  is  non-Lambertian  and  favoi  small  angles;  the  radiance  profiles  are  quite 
similar  to  each  other.  Figure  15  shows  a  composite  plot  of  all  the  transmission 
histograms  and  comparison  with  a  skewed-Lambertian  profile.  To  first  order,  it  can 
be  concluded  that  the  transmission  radiance  profile  is  independent  of  optical 
thickness  and  source  zenith  for  thick  clouds. 

The  skewed-Lambertian  radiance  profile  can  be  represented  bv  the  probability 
densitv  function 


p  ( 8  )  =  (  m  ■+■  1 

)cos  ’’t)  sin# 

!  9  ) 

or 

p ( u .  =  (m 

F 

a. 

+ 

■  10) 

where  m 
empirically 

=  I  tor  a  Lambertian  profile 
;  was  used  to  generate  'he  curve 

A  value  of  m  ~  1 
tor  Figure  15. 

! .  352 

■  determined 
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In  summary,  while  conventional  wisdom  holds  that  the  transmission  and 
reflection  functions  are  Lambertian,  thev  are  not.  Transmission  radiance  protiles 
are  always  skewed  to  low'er  angles  (i.e.,  they  are  more  forward  scattered  than 
Lambertian)  and  reflection  radiance  profiles  are  sensitive  to  cloud  optical  thickness 
and  source  zenith.  Only  for  reflection  of  diffuse  light  from  very  thick  clouds  <e.g., 
r  =  100)  was  a  Lambertian  radiance  protile  observed. 

3.3  Intensity  Reference  Method 

The  cloud  model  was  extended  to  include  a  surface  receiver  on  the  ground.  It  was 
pointed  out  earlier  that  to  actually  collect  photons  at  a  small  receiver  is  impractical 
and  alternative  methods  were  sought.  The  intensity  reference  method  (IRM), 
described  in  Section  2.2  and  Appendix  B  was  implemented.  Three  probability 
models  were  examined.  The  first  IRM  model  assumed  that  the  radiance  profile  at 
the  base  of  the  cloud  is  Lambertian.  In  that  case  tire  probability  of  a  photon  being 
scattered  to  the  receiver  is 


V  - 


(ID 


where  Aa>  is  the  solid  angle  defined  by  the  receiver  from  the  photon  exit  point  at 

the  cloud  base  and  u  is  the  direction  cosine.  The  solid  angle  is  written  as  a 
differential  as  a  reminder  that  the  analysis  is  valid  for  small  solid  angles.  For  larger 
solid  angles  the  probability  must  be  represented  a  an  integral,  viz. 

since  the  direction  cosine  cannot  be  regarded  as  a  constant  over  the  receiver  area. 

Table  4  shows  the  results  of  a  computation  with  a  Lambertian  IRM  model.  The 
results  are  highly  unsatisfactory  insofar  as  the  transmission  predicted  for  the  ideal 
receiver  is  24%  less  than  the  actual  value  (as  determined  from  the  actual  number  of 
photons  transmitted  through  the  cloud).  This  discrepancy  prompted  the  analysis  of 
cloud  radiance  profiles  discussed  earlier  in  Section  3.2.  It  was  concluded  that  the 
discrepancy  was  due  in  large  part  to  the  failure  of  the  Lambertian  radiance  profile  to 
adequately  predict  the  forward-scattered  nature  of  the  radiance  profile 

The  second  IRM  model  assumed  that  the  radiance  profile  at  the  base  of  the  cloud  is 
Msewed-Lambertian  (as  described  above)  In  that  case  the  probability  of  a  photon 
being  scattered  to  the  receiver  is 
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Table  4  shows  the  results  of  the  computation  with  the  skewed-Lambertian  IRM 
model.  The  results  show  considerable  improvement  with  the  error  being  brought 
down  to  lOT.  This  is  better  but  not  altogether  satisfying.  Recourse  was  made  to  a 
more  detailed  IRM  calculation 

In  a  formal  sense  the  IRM  should  be  based  on  the  actual  scattering  events  rather 
than  on  the  global  or  macroscopic  phenomena  such  as  radiance  profiles.  Thus,  the 
third  IRM  model  made  no  assumptions  about  the  radiance  profile.  Rather,  when 
the  photon  is  predicted  to  leave  the  cloud  (i.e.,  its  last  scattering  event),  the 
probability  of  the  photon  reaching  the  receiver  is  determined  by  the  scattering  phase 
function,  to  wit, 

p  =  P  {9  )d co  (14) 


where  P  (0)  is  the  normalized  phase  function.  The  Henvev-Greenstein  phase 
function  used  in  this  study  is  given  by 


P(0) 


g2 

(1+  g2~2gp) 


where  g  is  the  asymmetry  parameter. 

Table  4  shows  the  results  of  the  computation  with  the  Henvev-Greenstein  IRM 
model.  The  error  is  about  2%.  Clearly,  this  is  the  appropriate  IRM  model. 
Unfortunately,  it  is  somewhat  more  computation  intensive  than  the  other  methods 
because  of  the  scattering  geometry. 

The  Henvev-Greenstein  IRM  will  be  adopted  for  all  further  cloud  and  receiver 
simulation  in  this  report. 

3.4  RIMS  Simulation /Inferred  Optical  Thickness 

The  Remote  Irradiance  Measuring  Systems  (RIMS)  is  a  field  unit  for  measuring 
irradiance.  Optical  thickness  can  be  inferred  from  the  measured  irradiance.  The 
cloud  and  receiver  simulation  are  applicable  to  RIMS.  To  wit,  the  optical  thickness 
can  be  inferred  from  "measured"  (i.e.,  simulation)  photons  reaching  the  receiver. 
This  method  will  be  employed  to  determine  the  optical  thickness  in  model 
calculations  of  signal-to-noise  ratio. 
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3.5  Pulsewidth 

The  cubic  cloud  model  was  modified  to  record  the  temporal  behavior  of  the 
photons.  Specifically,  as  photons  move  through  the  cloud  and  to  the  receiver  the 
times  of  transit  me  accumulated  and  the  total  transit  time  statistics  are  collected. 
Two  types  of  statistics  are  developed.  One  is  the  simple  mean  and  standard 
deviation  while  the  other  is  a  histogram  of  transit  time. 

Figures  16-21  show  pulsewidth  histograms  for  the  ideal  and  actual  receivers  for  a 
variety  of  conditions.  Also  shown  for  comparison  is  the  Lee-Schroeder  pulsewidth 
model  (see,  for  example,  Lee,  et  al.  1986).  The  latter  model  consists  solely  of  a  3dB- 
pulsewidth  calculation,  but  a  modified  gamma  function  is  assumed  for  the  pulse 
shape.  Since  Figs.  16-21  represent  the  same  cloud  optical  and  physical  properties 

(r  =  20,  Z  =  H  =  1  km)  the  model  curve  is  the  same  in  all  the  figures.  These  results 
show  a  considerable  discrepancy  between  the  present  simulations  and  the  model 
The  present  results  show  a  dependence  of  the  pulsewidth  on  the  signal  spot  size, 
receiver  field  of  view,  and  receiver  ground  position  which  are  unaccounted  for  by 
the  model.  Generally  speaking,  the  model  tends  to  overestimate  the  pulsewidth 
and  the  results  get  worse  with  diminishing  receiver  field  of  view. 

Unfortunately,  the  model  was  never  properly  documented  and  so  its  origins  are 
unknowm.  Aiter  some  study  and  analysis  it  was  concluded  (i.e.,  guessed)  that  their 
model  was  based  on  a  point  source  of  photons  and  an  infinite  receiver  (i.e.,  the  time 
for  any  photons  reaching  the  surface  is  counted).  A  simulation  was  run  of  the 
presumed  model  described  above.  The  results  and  comparison  with  the  model  are 
shown  in  Figure  22.  These  results  are  in  very  good  agreement  with  each  other  and 
with  those  shown  in  Lee,  et  al.  (1986).  This  lends  credence  to  the  hypothesis  of  a 
point  source  and  a  distributed  receiver  and  may  explain  why  the  model  consistently 
overestimates  the  pulsewidth. 

3.6  Signal-to-Noise  Ratio  Simulation 

As  described  previously,  the  signal-to-noise  ratio  simulation  must  be  developed 
from  separate  simulations  for  the  signal  and  background.  A  suite  of  simulations 
was  established  for  comparison  with  the  link  model.  Tn  all,  twelve  simulations 
were  run  including: 

2  laser  zenith  angles, 

2  laser  spot  sizes. 

2  sun  angles, 

1  cloud  condition,  and 

2  receiver  locations. 
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labito  J  liuough  16  show  die  results  of  the  simulations.  The  analysis  ;s  limited  to  a 

single  cloud  condition  (r  =  20 ,  Z  =  H  =  1  km)  because  of  the  excessive  computing 
time  associated  with  the  simulations.  (Typically,  three  hours  are  required  on  a 
Compaq  386/20  witn  a  Weitek  coprocessor  to  get  a  decent  pulse  histogram.)  Also, 
the  receiver  field  of  view  half-angle  was  limited  to  45°.  Smaller  field  of  view 
translates  to  longer  computation  time  because  fewer  photons  land  in  the  field  of 

T.'i  jw 

Looking  at  Tables  5-16  the  following  is  observed:  The  user  supplies  the  solar  zenith, 
signal  nadir,  spot  radius  on  the  cloud,  cloud  top  area,  receiver  displacement  (along 
the  x-  or  major-axis),  the  cloud  thickness  and  height,  and  the  receiver  parameters 
used  in  the  simulation.  The  user  then  provides  the  numbers  of  photons  from  the 
results  of  the  background  and  signal  simulations  and  the  3dB  puisewidth  measured 
from  the  histograms  (Figures  16-21).  For  the  model,  the  spot  extension  must  be 
specified.  This  was  taken  as  approximately  unity  for  a  centered  receiver  and  one- 
half  for  a  receiver  placed  beneath  the  spot  edge.  The  optical  thickness  for  the 
downlink  calculation  is  inferred  from  the  background  simulation.  The  results  are 
shown  in  tine  lower  part  of  the  table.  The  3dB  puisewidth,  signal,  background,  and 
SNR  are  compared.  Also,  the  receiver  displacement  (simulation)  can  be  compared 
with  (half  of)  the  spot  major  axis  on  the  surface  (model). 

Generally  speaking  the  simulation  and  model  signal  predictions  are  in  g,Gcd 
agreement,  with  an  average  difference  of  about  1.1  dB.  There  are  only  two  different 
background  values,  and  they  agree  to  within  less  than  0.5  dB.  (These  results  would 
probably  have  been  better  if  a  larger  cloud  area  was  used,  but  that  would  have 
increased  computation  time).  The  signal-to-noise  ratio  for  the  simulation  was 
typically  5  dB  higher  than  that  of  the  model.  This  is  attributed  in  large  part  to  the 
discrepancv  of  the  simulation  and  model  puisewidth.  Figure  23  shows  a 
comparison  of  the  signal- to-noise  predictions  tor  the  simulation  and  model. 

4.C  CONCLUSIONS  AND  RECOMMENDATIONS 

A  Monte  Carlo  simulation  for  finite  parallelpiped  clouds  has  been  developed  and 
rested  for  consistency  with  stratus  (i.e.,  infinite)  cloud  models.  Radiance  profiles 
generated  by  the  simulation  were  shown  not  to  be  Lambertian  as  is  generally 
believed.  On  transmission,  the  radiance  profiles  appear  to  be  relatively  independent 
of  source  angle  and  optical  thickness  and  exhibit  a  skewed  radiance  protile  which 
*avors  smaller  angles.  On  reflection,  the  radiance  profiles  cannot  be  so  simply 
characterized.  Thev  are  probabiv  bidirectional  (i.e.,  dependent  upon  the  azimuth 
angle  as  well)  and  that  was  not  considered  in  the  present  report.  Only  in  the  case  of 
reflection  of  a  diffuse  ('Lambertian)  source  off  a  very  thick  cloud  was  a  Lambertian 
radiance  protile  observed. 
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Pulse  histograms  were  nerated  tor  a  number  of  source  and  receiver  combinations 
tor  a  given  cloud  condition.  The  puisewidth  and  shape  were  found  to  be  dependent 
unon  the  spot  size,  source  zenith,  and  receiver  field  of  view  and  surface  location 
beneath  the  cloud  in  addition  to  its  optical  and  physical  properties.  This  is  in  sharp 
contrast  with  the  model  currently  used  in  link  programs  which  depends  on.iv  upon 
the  optical  and  physical  properties.  Simulated  pulsewidths  were  of  the  order  of  one- 
half  of  the  predicted  values  for  the  actual  receiver. 

Signai-to-noise  ratios  generated  from  simulations  were  also  carried  out  and 
compared  with  model  results.  Generally,  the  background  values  agree  well  (less 
than  0.5  dB  difference)  and  the  signal  values  agree  well  (about  1  dB  difference).  But 
the  signai-to-noise  ratio  is  typically  about  5  dB  higher  in  the  simulation  than  in  the 
model,  due  in  large  par'  to  the  discrepancy  in  the  puisewidth. 

The  work  described  ;  t  this  report  shows  that  the  cubic  cloud  mode!  has  great 
rotential  for  expandi?  g  our  understanding  of  radiative  transfer  in  clouds  But  this 
:s  only  the  first  step.  This  model  is  a  valid  tool  for  studying  non-uniform  cloud 
optics  and  should  be  applied  to  the  dow’nlink,  uplink,  and  remote  sensing  problems. 
Also,  the  results  of  the  radiance  profile  and  puNewidth  indicate  that  the  spatial, 
angular  and  temporal  spreading  of  light  in  stratus  cloud  should  be  re-examined. 
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Table  1.  ResuL.  of  Cubic  Cloud  Simulation  (lxlxl  km  Cloud) 
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Table  2.  Results  of  Cubic  Cloud  Simulation  ox5xl  km  Cloud) 
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Table  3.  Results  of  Cubic  Cloud  Simulation  (10x10x1  km  Cloud) 
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Table  4.  Comparison  of  Cloud  Transmission  from  Intensity  Reference  Models 
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Figure  2.  Schematic  of  Intensity  Reference  Method 
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Figure  8.  Comparison  of  Simulation  and  Lambertian  Radiance  Profile 
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Figure  10.  Comparison  of  Simulation  and  Lambertian  Radiance  Profile 
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Figure  14.  Comparison  of  Simulation  and  Lambertian  Radiance  Profile 
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APPENDIX  A.  MONTE  CARLO  SIMULATION  OF  RADIATIVE  TRANSFER 
1.0  INTRODUCTION 

The  purpose  of  this  Appendix  is  to  give  a  brief  introduction  to  and  description  of 
Monte  Carlo  simulation  methods  for  radiative  transfer  problems.  The  Appendix  is 
intended  as  an  adjunct  to  the  Monte  Carlo  codes  developed  for  the  Naval  Ocean 
Systems  Center;  its  function  is  to  define  the  underlying  basis  of  the  codes  and  does 
not  purport  to  be  a  primer  on  Monte  Carlo  simulations.  The  reader  is  referred  to 
Kalos  and  Whitlock  (1986)  for  a  general  discussion  of  Monte  Carlo  methods  and  to 
Lenoble  (1985)  for  a  discussion  of  its  application  to  radiative  transfer. 

In  a  Monte  Carlo  computation  one  photon  at  a  time  is  followed  along  its  three- 
dimensional  path  through  a  scattering  medium.  Its  fate  is  determined  by  suitable 
probabilitv  distributions  for  mean  free  path,  absorption,  scattering  angle,  wall 
absorption  and  reflection,  and  so  on.  The  photon  is  followed  until  it  is  absorbed, 
detected,  or  leaves  the  field  of  interest.  A  sufficient  number  of  photons  are  followed 
until  a  picture  of  the  system  emerges. 

The  Monte  Carlo  method  has  some  advantages  over  other  computational  methods 
for  radiative  transfer,  namely 

•  any  phase  function  can  be  used 

•  can  include  polarization  (with  a  two  times  penaltv  in  computation  time) 

•  several  detectors  may  be  included  (small  penalty  in  computation  time) 

•  can  divide  medium  into  both  vertical  and  horizontal  layers  of  different 
optical  properties  (small  penalty  in  computation  time) 

•  complex  geometries  can  be  used  (small  penalty  in  computation  time) 

Naturally,  there  are  disadvantages  as  well  .  .  . 

•  not  suitable  if  high  accuracy  is  desired  (doul  mg  the  accuracy  requires 
quadrupling  the  computer  time) 

•  always  get  bin  averaged  radiances  rather  than  point  values  (smaller  bins 
mean  increased  computation  time  for  a  given  accuracy) 

•  not  practical  for  opticallv  thick  media  (say,  r  >  100) 


A  technical  discussion  of  the  Monte  Carlo  method  for  radiative  transfer  is  presented 
in  Section  2.  The  phvsics  is  described  and  the  geometry  of  scattering  is  presented. 
Next,  the  statistical  relationships  required  are  reviewed.  Then  the  probability 
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distributions  for  mean  free  path  and  scattering  angles  are  discussed.  Finally,  some 
miscellaneous  subjects  such  as  absorption  and  reflection  are  discussed. 

Section  3  contains  a  discussion  of  the  computational  aspects  of  Monte  Carle 
simulation.  Program  design  philosophy  and  layout  are  discussed  followed  bv  a 
collection  of  tricks  the  author  has  found  useful  for  reducing  the  computation  time. 

References  and  List  of  Figures  are  gi”en  in  Sections  4  and  5,  respectively. 

2.0  TECHNICAL  DISCUSSION 

Figure  A-l  shows  a  schematic  of  the  scattering  process.  Photons  are  introduced  into 
the  scattering  medium  from  a  source  with  a  specified  mean  free  path  and 
orientation.  Upon  encountering  a  scattering  particle  the  photon  is  scattered  off  with 
a  new  mean  free  path  and  direction.  When  the  photon  is  scattered,  its  departure 
direction  is  determined  relative  to  its  arrival  direction.  In  order  to  properly  follow 
the  photon  its  position  in  absolute  coordinates  must  be  determined. 

2.1  Geometrical  Considerations 

Figure  A-2  shows  a  sketch  of  the  geometry  of  scattering.  The  primed  coordinate 
system  is  relative  to  the  arrival  direction  (i.e.,  z'  aligns  with  the  direction  of  the 
arriving  photon).  The  unprimed  coordinate  system  is  aligned  with  the  absolute 
coordinates. 


Consider  a  scattering  event  with  a  unit  mean  free  path  scattered  at  a  zenith  angle, 

O'  and  azimuth  angle,  a'  relative  to  the  arrival  direction.  Then  the  following 
transformation  of  coordinates  applies 

/  / 

'a  *  * 
x'  ~  sincfc'cos#' 

f 

y'  =  sina'sin#  (A-l) 

t 

z'  =cosO 


x  —  —  u'  sin  a'  +  x"  cos  a' 

y  -  x " s i n a '  +  y'cos# 

/  / 

z  -  -  r '  sint?  -r-  z  '  cos# 


(A- 
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and 


9  =  cos  z 


ar  =  tan"’ 


(A-4) 


Equations  (A-l)  through  (A-4)  must  be  computed  for  each  scattering  event.  This  can 
literally  be  in  the  tens  of  millions  for  a  modest  simulation.  Unfortunately, 
computation  times  can  be  quite  large  because  of  the  trigonometric  functions.  After 
much  experimentation  with  alternative  computation  schemes  (including  fast  look¬ 
up  tables)  it  was  found  that  the  entire  calculation  could  be  done  much  faster  in 
terms  of  direction  cosines,  and  that  very  few  trigonometric  calculations  are  required 
at  all.  In  fact,  as  will  be  seen  later  on,  even  the  scattering  angle  zenith  can  be 
randomly  drawn  in  terms  of  its  direction  cosine.  (The  azimuth  is  still  drawn  as  an 
angle,  however.) 


Referring  to  Figure  A-3,  the  direction  cosines  of  the  translation  vector  ix  from  the 
point  of  scattering,  x,  to  the  next  collision  J  are  given  by 


where 


1  "  /  "l  ' 

,  lS~l2l2 

V’-c 

(AS) 

.  //  .  //,  / 

,  l2C  +  l:l2 

(A -b) 

h  7'- !? 

/  rr  2  t  ff  / 

1} =  "  V  1  “  1  3  l:+  1  2l  2 

(A  -7) 

( A-8) 

and  I,  are  the  direction  cosines  of  the  scattered  photon  in  absolute  coordinates,  /, 

// 

are  those  in  coordinates  relative  to  the  arriving  photon,  and  l  .  are  those  of  the 
arriving  photon  in  absolute  coordinates.  Notice  that  there  is  no  explicit  dependence 
on  anv  angles. 


For  a  mean  free  path  of  length  d  the  final  position  x.  is  determined  from 


Optical  Transmission  Through  Clouds 

Appendix  A.  Monte  Carlo  Simulation 

Page  A-4 


It  has  been  shown  that  the  photon  position  can  be  tracked  bv  calculating  its  position 
from  its  previous  position  and  translation  vector.  The  translation  vector  is  obtained 
from  a  randomly  drawn  mean  free  path,  direction  cosine  of  zenith,  and  azimuth 
angle. 

2.2  Statistical  Considerations 


Referring  to  Figure  A-4,  p(x)  is  called  a  probability  density  function  and  its  integral 
hr)  is  called  a  probability  distribution  function.  The  definition  of  the  probabilitv 
densitv  function  requires  that 


p(x)  >0 
j  p  (x)dx  =  1 


(A- 10) 


The  probability  distribution  function  is  a  cumulative  probability,  i.e.. 


P  (X  <,  x  )  =  f(x  )  =  j  p  ( t  )d t  ( a-  1 1 ) 

The  distribution  function  satisfies  these  three  conditions: 

1.  \imf(x)  =  Q;  lim/(r)  =  l 

x  -o  1  — 

2.  f  <  x  )  >  0;  /  (  y  )  >  f  (x  )  if  y  >  x 

3.  fix  )  is  continuous 

Finally,  to  get  a  random  variable  x  with  a  distribution  function  *'x),  choose  a 
random  number,  RX,  uniform  in  [0,1]  and  get  x  from  the  inverse  function,  i.e., 

x  =  r'(RN)  ( A- 1 2) 


always  worthwhile  to  test  that  an  algorithm  indeed  samples 


hr).  Several 

methods  to  check  an  algorithm  can  be  used.  The  simplest  consists  of  generating 
random  variables  and  sorting  the  results  of  Eq.  1A-I2)  into  bins  within  the  range  of 
the  random  variable.  This  can  then  be  compared  with  the  probabilitv  density 
functmm  Several  euuu.jT.-"  w«.  oe  ^^.r.  Idter  ,n  connection  with  distribution 
functions  of  interest  m  radiative  transfer. 
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2.3  Mean  Free  Path 

Simulation  of  the  distance  traveled  between  collisions  is  calculated  directb'  since  the 
fraction  of  radiation  transmitted  through  a  given  distance  is  also  the  probability  of  a 
photon  traveling  the  same  distance.  Thus, 


f 

V  ( x  )  =  exp  ( -r )  =  exp 

v 


J  P*s 

0  J 


1  A- 1 3) 


f  ( a  )  =  1  -  exp  ( —x ) 


CA-14) 


^pds  =  r  =  -  1  n(  1  -  RN)  (A-15) 

o 


If  the  scattering  coefficient  is  constant  then  the  mean  free  path  is  given  simply  by 


s  =--jln(RN) 


(A-16) 


Figure  A-5  shows  a  comparison  of  the  statistical  model,  Eq.  (A-16)  with  the 
probabilitv  density  function.  Eq.  (A-13)  (for  a  constant  scattering  coefficient).  This 
demonstrates  that  applying  Eq.  (A-16)  to  random  numbers  uniform  in  [0,1]  does 
indeed  replicate  the  expected  probability  density  function. 

More  generally,  Eq.  (A-15)  must  be  iterated  to  get  the  mean  free  path.  s.  For  a 
simple  two-layer  model,  however,  analytic  results  can  also  be  obtained.  The 
important  point  there  is  to  adjust  the  mean  free  path  if  the  photon  goes  from  a 
region  of  one  optical  density  to  the  other.  Figure  A-6  shows  a  sketch  of  the  two-laver 

model. 

Equation  :  A-15)  suggests  that  a  medium  ot  continuously  variable  optical  properties 
might  be  treated  simply  bv  transforming  from  the  physical  plane  to  that  of  the 

ipticai  thickness,  r.  This  has  never  been  tried  by  the  author. 


Scattering  Functions 


chan 


scatter 


different  scattering  functions  have  been  used  in  Monte  Carlo  Simula 
ve  transport  for  different  physical  problems.  Figure  A -7  shows  a  skete 
;ng  geometry.  The  brief  list  below  is  limited  to  those  with  whicn  the  a 


tions  of 
h  of  the 
uthor  is 


familiar: 
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•  uniform  distribution  (usually  applicable  to  the  azimuth  angle  of 
scattering) 

•  isotropic  (equal  distribution  of  photons  per  unit  solid  angle;  usually 
application  to  photon  emission  on  quantum  state  transitions) 

•  Lambertian  (equal  distribution  of  photons  per  unit  solid  angle  per  unit 
area;  usually  applicable  to  reflection  from  surfaces) 

•  Mie  scattering  (complex  scattering  pattern  caused  by  the  interaction  of 
electromagnetic  waves  with  molecular  dipoles  characterized  by  strong 
forward  scattering,  usually  applicable  to  scattering  in  aerosols  and  clouds). 
Figure  A-8  shows  a  schematic  of  Mie  scattering 

•  Henyev-Greenstein  (an  analytical  model  with  characteristics  of  Mie 
scattering) 

•  Irvine-Henvey-Greenstein  (an  extension  of  the  above  model  with 
improved  backward  scattering) 

•  Empirical  (phase  function  and/or  distribution  are  fit  to  experimental  data; 
these  have  been  used  for  clouds  and  sea  water,  for  example) 

•  Skewed-Lambertian  (cosine-to-the-n  behavior) 

Other  possibilities  abound  as  well.  King  and  Harshvardhan  (198b)  and  others  have 
used  phase  functions  derived  from  Mie  scattering  calculations  for  a  cloud  with  a 
particle  size  distribution.  Waldman  (1988)  has  used  curve  fits  to  experimental  data 
of  Petzold  (1972)  to  get  the  scattering  function  in  sea  water.  Some  specific  examples 
are  considered  in  detail  below. 

A.  Uniform  Distribution 

An  example  of  uniform  distribution  is  the  azimuth  angle  upon  scattering.  The 
angle  is  selected  randomly  from 

0=2 n RN 


where  RN  is  uniform  in  (0,1]. 

8.  Isotropic  Distribution 

Isotropic  distribution  applies,  for  example,  to  spontaneous  emission  of  a  photon 
during  a  change  in  quantum  state  such  as  fluorescence.  Statistically,  the  intensity,  or 
number  of  photons  per  steradian  is  constant,  therefore 
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p,(9)  =  sin  9 

f[(Q)  =  1  -  cost?  (A-17) 

p  =  RN 

(Notice  that  RN  is  functionally  equivalent  to  1-RN  since  RN  is  uniform  in  [0,1].) 

C.  Lambertian  Distribution 

This  is  applicable  to  diffuse  reflection  from  a  surface.  Statistically,  the  radiance  or 
number  of  photons  per  unit  area,  per  steradian  is  constant,  therefore 

pL(0)  =  2sin 9  cos 9 

fL(9  )  =  1  -  cos  70  (A-18) 

p  =  a/rN 

Figure  A-9  shows  a  comparison  of  the  statistical  model  and  the  probability  densitv 
function. 

D.  Henyev-Greenstein  Distribution 

This  is  a  commonly  used  model  for  Mie  scattering.  Although  it  lacks  the  details  of 
the  Mie  scattering  it  has  the  general  characteristics  and  often  proves  to  be  an 
adequate  model  for  cloud  and  aerosol  studies.  Its  principal  advantage  is  its  analytic 
form,  viz.. 


Phc 


(60  = 


(1  -  g2)sin6 


2(1  +  g2  -  2 g  cos 9  ) 


3/2 


,(«,)  =  I  !uf_Y_! _  i 

M  \  ^  %  -J  1  +  gz  ~  2 g  cos  9  J 


(A- 19? 


1  +  g2  - 


i-  r 


p  = 


,  1  +  g  -  2  g  R  N  ) 


Figure  A- 10  shows  a  comparison  of  the  statistical  model  and  the  probability  densitv 
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E.  Irvine-Henyey-Greenstein  Distribution 

The  principal  shortcoming  of  the  Henvey-Greenstein  phase  function  is  its  weak 
backward  scattering  lobe.  The  Irvine  model  attempts  to  correct  this  bv  summing  two 
Henvey-Greenstein  phase  functions  with  positive  and  negative  asymmetries, 
respectively.  Thus, 


P:HG(0)  =  apHG(9;g,)  +(1-  a)pHG(d;g2) 

fine (d  >  =  af»c  (6  -Z  i)  +  (1  -  a  )fHC  (6 ;  g2)  '  A'20/ 

where  a  is  the  fraction  of  forward  scattered  photons.  Generally,  g,  >  0  and  g,  <  0 . 
The  direction  cosine  must  be  obtained  from  the  distribution  function  by  iteration. 
Figure  A- 11  shows  a  comparison  of  the  statistical  model  and  the  probability  density 
function. 

F  Sea  Water  Scattering  Distribution 

The  Henvey-Greenstein  phase  function  is  not  generally  adequate  to  model  scattering 
of  light  in  sea  water  because  a  high  absorption  coefficient  puts  a  stronger  emphasis 
on  the  forward  and  backward  scattering  lobes.  Statistical  models  for  the  sea  water 
scattering  distribution  function  can  be  developed  from  empirical  data,  such  as  that 
presented  bv  Petzold  (1972). 

Figure  A-12  shows  a  comparison  of  such  a  statistical  model  (curve  fit  to  data)  and 
the  probability  density  function  (data  on  the  volume  scattering  function  in  sea 
water). 

G.  Skewed-Lambertian 

There  are  no  specific  applications  for  this  distribution.  Frequently,  however,  one 
encounters  radiance  profiles  which  are  described  as  having  "cosine-to-the-ri 
behavior,  i.e.. 


p  ..  id  )  =  (  n  +  I) cos  " 9  sin d 

r  'ft, 

q.  [0  )  =  1  -  cos  ’’*  0 
a  -  R  N"-" 


mows  a  comparison  of  the  statistical  model  and  the  probability  density 
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2.5  Surface  Reflection 

Surface  reflections  mav  be  diffuse  or  specular.  One  special  case  of  diffuse  reflection 
is  Lambertian,  in  which  the  reflected  radiance  is  the  same  in  all  directions.  Figure 
A- 12  shows  a  sketch  of  the  reflected  irradiance  for  various  types  of  surface  reflection. 

A  commonlv  encountered  problem  in  radiative  transfer  is  the  specular  reflection  of 
photons  off  curved  surfaces.  Here  the  vector  form  of  Fermat's  principle  is  required 
to  determine  the  direction  of  the  reflected  photon, 

r=  f  -  2(n  ■  r~)n  CA-22) 

where  f  and  r"  are  the  incident  and  reflected  directions,  respectively,  and  n  is 
the  unit  normal  to  the  surface  (see  Figure  A-12). 

When  the  walls  have  discontinuous  first  derivatives  (i.e.,  corners,  as  in  the 
boundary  of  a  wall  and  floor)  then  special  attention  must  be  paid  to  the  photon 
reflection.  As  unlikely  as  such  an  encounter  may  be  in  a  simulation,  experience  has 
shown  that  it  can  and  will  happen,  often  with  dire  consequences  in  a  program 
which  did  not  anticipate  it. 

2.6  Absorption 

Absorption  can  occur  within  the  scattering  medium  or  at  a  surface.  In  either  case 
the  probability  of  the  photon  being  absorbed  can  be  computed  and  the  photon  path 
can  be  continued,  or  not,  depending  on  the  outcome. 

This  can  be  verv  inefficient  from  a  computational  point  of  view,  however.  Rather, 
what  is  recommended  is  that  each  photon  is  a  'signed  a  weight,  or  probability,  which 
is  diminished  upon  each  absorption.  The  decrease  in  photon  weight  is  just  the 
probability  chat  it  was  absorbed.  In  this  way  the  computation  can  proceed  without 
starting  a  new  photon  each  time  one  is  absorbed. 

2.7  Gaussian  Radiance  Profile 

This  special  topic  is  of  relevance  to  photon  detection  under  water  where  the 
radiance  profile  is  tvpicallv  Gaussian  (see,  for  example,  Jerlov,  W76>.  There  are 
manv  wavs  to  generate  random  variables  which  replicate  Gaussian  probability 
density  functions,  but  these  density  functions  are  not  the  same  as  those  associated 
with  a  Gaussian  radiance  profile  The  difference  is  that  the  radiance  protile 
ton  tarns  an  additional  cosine  term  '  for  the  area  effect)  and  a  sine  term  'tor  the  solid 
angle  effect  >. 
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In  a  formal  sense  the  probability  density  function  for  a  Gaussian  tadiance  profile  for 
a  specitic  receiver  is  given  by 


|  exp  (-(O  cosd  sin 9  d  6 

r  <#)=-£ - ( A-23) 

I  exp^-(lr)  ycose?  sin 0d6 

"  r\ 


where  0.,r  is  the  receiver  field  of  view  half-angle,  it  has  been  assumed  that  the 
receiver  has  a  cosine  response  (otherwise  substitute  the  response  function  for 

COS d  )■ 

The  probability  of  integrating  this  in  closed  form,  then  determining  the  distribution 
function  and  its  inverse,  is  zero.  Waldman  (1987)  has  found  a  suitable 
approximation  for  the  exponential  term  which  permits  integration  of  Eq.  (22). 
Nevertheless,  the  distribution  function  is  cumbersome  and  the  direction  cosine  can 
only  be  obtained  by  iteration.  This  is  somewhat  time  consuming  and  a  table  look-up 
is  prefer. ed.  Figure  A-13  shows  a  comparison  of  the  statistical  model  and  the 
probability  density  function  for  a  receiver  with  a  90°  field  of  view  half-angle  (fiat 
plate  receiver). 

2.S  Monte  Carlo  Simula  bon  Results 

Several  kinds  of  results  can  be  gleaned  from  Monte  Carlo  simulations.  These  can  be 
broadly  characterized  as  global,  statistical,  and  specific  As  an  example,  consider  light 
iransmission  from  a  point  source  through  a  cloud-  Global  results  consist  of  the  total 
numoer  of  photons  reflected  and  transmitted.  Statistical  information  might  consist 
of  .he  mean  and  rms  values  for  exit  radius,  angle,  and  time.  The  mean  and  rms  are 
raicuiated  as  follows: 


V 


A -24) 


■j'  - 


Y 


A-25) 
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y  —  -  v  y  v  y  >  -A-2bi 

where  ?<  are  the  functions  being  evaluated  (e  g...  exit  radius,  angle,  or  time)  and  a/, 
is  a  weigh  ting  function.  The  weighting  function  applies  in  such  cases  where  the 
photons  may  arrive  at  the  exit  plane  with  different  weights  (e.g.,  due  to  absorption,' 
These  statistics  can  be  computing  on  a  running  basis  so  that  it  is  not  necessary  to 
store  ad  or  the  individual  results  (which  would  be  impractical,  if  not  impossible,  on 
most  computers). 

Specific  results  refer  to  such  information  as  the  radiance  profile,  surface  energy 
distribution,  and  pulse  temporal  distribution.  The  nature  of  the  Monte  Carlo 
simulation  is  that  these  cannot  be  determined  exactly,  but  rather  only  as  bin 
averaged  quantities.  To  collect  this  information  requires  separating  the  results  te.g  , 
exit  angle.?  into  discrete  bins  and  accumulating  the  weighted  photons  that  arrive  in 
each  bin.  These  data  can  then  be  plotted  as  histograms  to  give  an  approximate 
picture  of  the  photon  behavior.  Clearly,  this  contains  the  most  information  of  ail 
the  results  but  can  be  particularly  time  consuming,  depending  upon  the  bin 


C  OM  PUT  .ATI OX  A  L  A  S  P ECTS 


Monte  Ca 
solutions 

meet  ins  a 


no  simulations  oniv  exist  because  of  computers,  there  are  no  analytical 
Program  design  and  Lvout  will  dictate  how  useful  the  program  Is  in 
variety  of  needs.  Program  execution  speed  is  vital  in  determining  if  the 


3 1  <3li. 


■am  Design  Philosophy  and  Layout 


vpehv  designed  Monte  Carlo  simulation  program  should  be  readily  adaptable  to 
number  of  -adiative  transfer  problems  within  a  broad  class  of  such  problems, 
-ome  event  this  means  anticipating  which  parameters  are  fixed  and  which  are 
a  Pie  Also  the  program  should  be  highly  modularized  so  that  anv  part  can  be 
:;!v  identified  and  altered  an  needed  Seif-documenting  programs  can  be  verv 

■  to  maintain  and  modify  This  can  be  achieved  bv  assigning  appropriate  ’'aria pie 
tes  that  are  either  descriptive  or  conform  to  the  standard  not  a  non  of  t  he  held, 
j,  identifying  ail  the  constants,  parameters,  and  variables.,  and  their  units,  along 

■  a  modest  amount  of  comments  ,n  the  .ode  is  indispensable 
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take  advantage  of  faster  machines  and  compilers  as  they  become  available.  During 
program  development  it  is  desirable  to  test  modules  independently  and  also  to  build 
in  debugging  statements  which  can  be  turned  on  at  will  .  (Actually,  with  modem 
compilers  and  debuggers  this  is  probably  not  necessary.)  Finally,  the  programming 
should  be  kept  simple.  There  are  times  when  this  will  be  in  direct  conflict  with 
making  the  program  as  fast  as  possible.  Of  course,  make  the  program  run  faster,  but 
document  the  simple  procedure  and  the  origin  of  the  more  arcane  one  in  the  code 
itself. 


3.2  Program  Execution  Speed 

in  Monte  Carlo  simulations  program  execution  speed  mav  mean  the  difference 
between  getting  an  answer  or  not.  Methods  of  increasing  program  speed  fall  into 
two  categories,  general  and  ad  hoc.  General  methods  consist  of  programming 
practices  which  are  germane  to  ail  simulations.  For  example,  it  has  already  been 
mentioned  that  trigonometric  calculations  are  slow  compared  to  those  for  direction 
cosines.  This  applies  to  scattering  statistics  as  well,  where  the  direction  cosines  can 

be  found  directly.  Along  the  same  lines,  yf \  -  /r  is  much  faster  than  sine  function; 

just  beware  of  angles  between  n  and  2 re  .  Frequently,  look-up  tables  can  be  faster 
than  computation  and  can  contain  enough  entries  so  that  interpolation  is  not 
required.  Using  pointers  to  the  table  (available  in  the  Pascal  and  C  languages)  is 
even  raster  Avoid  loops,  particularly  for  updating  direction  cosines  in  position. 
There  is  a  large  computational  overhead  involved  with  a  loop.  Finally,  use  the 
fastest  floating  point  representation  available  on  the  target  machine  and  compiler; 
there  are  no  hare  and  fast  rules  on  this.  On  some  machines /compilers  reals  are 
raster  than  doubles  because  the  data  fetches  are  quicker,  while  on  others,  thev  are 
slower  because  the  numerical  coprocessor  must  convert  all  the  reals  to  doubles 
anyway’  Of  course,  these  considerations  mav  be  overridden  bv  questions  of 
r.  jmericai  precision  and  range. 


When  dealing  with  specific  problems  some  ad  hoc  tricks  mav  be  found  that  can 
reduce  computation  time  significantly.  Without  any  elaboration,  these  mav  consist 
:f  -elective  generation  of  photons  (i.e.,  ignoring  ones  which  initially  go  off  in  the 
wrong  direction),  discarding  photons  which  stray  too  far,  calculating  several  cases  at 
on ce  if  possible  (e.g.,  several  receiver  locations,  different  wavelength  photons,  etc.), 
empdfv  complex  media  (e.g.,  in  an  atmospheric  transmission  problem  treat  a  cloud 
as  i  diffuse  reflector  rather  than  simulating  it),  exaggerate  receiver  si/e  to  collect 
more  rho‘ons,  Meat  on  receiver  shape  to  simplify  detection  calculation,  and  so  on. 


:g  up  the  calculations  computes  the 
•  er  than  yjnt  collecting  those  photons 
a  ole  and  ,;ome  cases  and  not  apphoat 
rer.ee  method  and  it  is  discussed  s»  pa 
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X 


To  get  a  random  variable  x  with  the 
distribution  function  f(x),  choose  a  random 
variable,  RN,  uniform  in  [0,  1]  and  get  x  from 
the  inverse  function,  i.e.. 


x  =  f  1  (RN) 


Probability  Density  Function 
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Exponential  (Poisson)  Distribution 


Length, Mear 


?e  ^ath 


isolation  73.  Equation 
iC  000  Random  Trials 


Figure  A-5  Comparison  of  Simulation  and  liquation: 
Exponential  (Poisson)  Distribution 


Simulation  of  the  distance  traveled  between  collisions  is  calculated  directly 
since  the  fraction  of  radiation  transmitted  through  a  given  distance  is  also  the 
probability  of  a  photon  traveling  the  same  distance 
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TWO-LAYER  MODEL 
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Figure  A -7.  Scattering  Geometry 
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Note: 

Complex  scattering  pattern  is  caused  by  the  interaction  of 
e'ectromagnetic  waves  with  molecular  dipoles  characterized  by 
strong  forward  scattering. 


Figure  AS.  Schematic  of  Mie  Scattering 
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Figure  A- 10.  Comparison  of  Simulation  and  Lquatu-n. 
I  lenvev  -Groenstein  Distribution 
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Skewed-Lambert ian  Distribution 


Simulation  v 3 ,  Equation 
130,000  Random  Trials 
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A-I3.  Comparison  of  Simulation  and  Equation: 
Skewec  '.ambertian  Distribution 
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Figure  A-14.  Surface  Reflection 
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Figure  A-I5.  Comparison  of  Simulation  and  Equation: 
Gaussian  Radiance  Distribution 
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APPENDIX  B.  INTENSITY  REFERENCE  METHOD  FOR  RADIATIVE  TRANSFER 
1.0  INTRODUCTION 

The  purpose  of  this  Appendix  is  to  give  a  brief  description  of  the  intensity  reference 
method.  The  method  has  been  described  by  Meier,  Lee,  and  Anderson  (1978)  and 
will  be  reviewed  briefly  here  along  with  a  description  of  a  model  extension  for  cases 
where  the  photon  is  near  the  receiver.  Basically,  there  are  three  distinct  domains  of 
interest  which  are  referred  to  as  far-field,  mid-field,  and  near-field  depending  upon 
the  distance  of  the  photon  to  the  detector  The  technical  discussion  contains  a  Liief 
description  of  each  of  these. 

In  the  intensity  reference  method  the  probabilitv  of  reaching  the  detector  is 
computed  and  stored  at  each  photon  scattenng.  The  accumulated  probabilities  are 
then  proportional  to  the  intensity,  as  described  below-.  Photons  continue  to  follow 
their  natural  histories  as  determined  by  the  specified  physical  processes.  The 
advantage  of  this  scheme  is  that,  in  contrast  to  the  low  likelihood  of  a  photon 
actually  hitting  a  receiver,  there  will  always  be  some  probability  of  it  actually- 
occurring.  Thus,  summing  the  probabilities  greatly  improves  the  statistics  since 
each  scattering  event  contributes  to  the  intensity.  In  addition,  the  probability-  is 
computed  exactly  from  the  scattering  point  to  the  receiver;  no  other  approximations 
are  introduced.  Moreover,  the  technique  lends  itself  to  vignetted  receivers,  i.e., 
receivers  with  a  limited  field  of  view  (w-rtich  lowers  the  likelihood  of  a  photon 
hitting  the  detector  even  further). 

2.0  TECHNICAL  DISCUSSION 

The  intensity  reference  method  considers  the  probability  of  the  photon  reaching  a 
detector  as  a  consequence  of  a  scattering  event.  Figure  B-l  -.hows  a  schematic 
drawing  of  the  method.  The  probability  of  a  photon  being  scattered  in  the  direction 
of  the  detector  is 

p  (6L)  =  [P  (6,  )da> 

(B-D 


where  the  integration  is  bounded  by  the  solid  angle  from  the  scattering  event  to  the 
receiver  area.  However,  along  the  path  to  the  detector  the  photon  may  be  absorbed 
or  scattered  out  of  the  path.  Thus,  the  probability  of  the  photon  actually  reaching 
the  receiver  is  diminished.  In  addition,  the  photon  "weight'  upon  scattering  is 
accounted  for  so  that  the  probability-  is  finally  given  bv 
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p  =  \wP  (Oje  "  4  d  co 

C  (B-2) 

where  w  is  the  initial  photon  weight,  a  and  3  are  the  absorption  and  scattering 
coefficients,  respectively,  and  d  is  range  from  the  scatterer  to  the  differential  solid 
angle.  Equation  (B-2)  is  a  two-dimensional  integration  and  is  numerically 
cumbersome  . 


2.1  Far-field  Approximation 

The  far-field  approximation  is  frequently  employed  in  radiative  transfer  problems. 
It  is  applicable  when  the  detector  area  is  small  compared  with  the  square  of  the 
range;  in  other  words,  for  small  solid  angles.  In  that  case,  Eq.  (B-2)  becomes 

p  =  wP(e,)s-<a^)t/CM;  dw 

where 

A  r  cos  £  Ar  cos  3£ 

______ 


(B-3) 

(B-4) 


and  9,  is  the  scattering  angle  given  by 

cos#,  =  I  l2  (B-5) 

where  h  is  the  direction  cosine  (unit  vector)  of  the  scattering  point,  x.  and  L  is  the 

direction  cosine  of  the  photon  arriving  at  x.  .  This  is  the  most  common  application 
of  the  intensity  reference  method.  It's  particularly  useful  in  atmospheric  problems 
where  the  distances  are  large  and  the  receivers  are  small. 

2.2  Mid-field  Calculation 

The  far-field  approximation  breaks  down  when  the  photon  is  sufficiently  close  to 
the  receiver  because  the  solid  angle  can  no  longer  be  approximated  by  Eq.  t B-4) . 
Thus,  recourse  must  be  made  to  Eq.  (B-!)  For  a  circular  receiver  Eq.  <B-1)  can  be 
written  as 


V 


wPlOje1"^ 


cos  C 


r  dr  d9 
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d 


(B-6) 
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where  R  is  the  receiver  radius.  The  distance.,  d.  the  scattering  angle,  Bi .  and  the 
zenith  angle,  C.  are  ail  functions  of  the  independent  variables  r  and  0 
Comparisons  of  computations  with  Eqs.  (B-3)  and  t B-6 )  show  that  Eq.  (B-3)  is  suitable 
if  the  range  :s  about  ten  times  the  receiver  radius.  Generally  speaking,  Eq.  (B-6)  is 
not  suitable  for  inclusion  in  Monte  Carlo  simulations  but  should  not  be  ruled  out 
altogether. 

1.3  Near -field  Approximation 

As  the  photons  get  closer  to  the  receiver  the  accuracy  required  of  the  integration 
increases,  thus  slowing  it  down.  This  occurs  at  a  height  above  the  -eceioer  of  about 


m  f-h  Af  ‘  hp 


one- tenth  or  me  melius.  Jn  that  close  r &ng6  the  solid  dngie  is  very  nearly  2z 
half  of  the  total  solid  angle).  Thus,  the  integral  can  be  approximated  as  follows: 


where  absorption  and  scattering  have  Peer,  neglected  (the  intensity  reference 
method  wouidn  t  make  seme  if  the  medium  was  so  optically  dense  at  that  distance 

‘o  the  receiver;,  k  is  the  direction  cosine  of  the  arriving  photon,  and  t  and  r  ere 
the  forward  and  backward  scattering  functions,  viz., 

f  \P-.vdn 


b  ~  2  ,T  1  P  i  a  a, 


B-~  is  clearlv  correc*  when  the  photon  is  iieaded 

•o  might  down  (l  =  -1,  p  ■-  :uf),  and  straight  act  ,-ss  : 
:Piv  reasonable  tor  anv  value  oi  direction  cosi  :e. 


CONCLUSIONS 


The  .r.ter.sitv  reference  met  hod  is  a  poweriu.  tool  for  Monte  (.  irio 
■mould  oe  used  where  :t  is  applicable. 
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Figure  B-I.  Schematic  of  Intensity  Reference  Method 
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APPENDIX  C.  MODIFIED  GAMMA  FUNCTION  PULSE 


The  modified  function  is  frequently  used  as  a  model  for  the  temporal  pulse 
spreading  of  scattered  radiation.  The  purpose  of  this  Appendix  is  to  collect  the 
appropriate  relations  for  the  pulse  time  to  peak,  mean,  rms,  and  3dB  pulsewidth. 
The  pulse  is  expressed  as  a  probability  density  function 

pit)  -  k2t  e~kt  (C-l) 

•where  k  is  the  inverse  characteristic  time  The  distribution  function  is  then 


f  ( t )  =  l  pit)  dt  =  I  -  ( 1  +  k  t  )e  ~t! 

0 

so  that  hQi  a  0  and  /( )  =  1.  The  time  to  peak  is  given  by 


The  mean  time  is  given  by 


dp  it  } 
dt 


*0:  O..  =  f 


le  rms  time  is  given  by 


f  =  Jtpit )dt  =  77 


U.  -  \(t  ~  t  )  pit  'd 


i  t  — 


\f'2 


U.  =  f(f  )  p  'i  '■it  =  ~~ 


«C-2) 


(C-3) 

(C-4) 


(C- 


( C  -6 ) 


The  3dB  ruise  width  is  given  hv  the  difference  in  times  where  pit !  -  piyii/2;  this 
o  ust  be  solved  bv  iteration,  there  results 


2.44638^03925705 
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